﻿#include "PSSDoseAveragedLET.hh"

#include <cmath>

#include "G4SDManager.hh"
#include "G4StepPoint.hh"
#include "G4SystemOfUnits.hh"
#include "G4TouchableHistory.hh"
#include "G4UnitsTable.hh"
#include "G4VPhysicalVolume.hh"
#include "G4VSensitiveDetector.hh"
#include "G4VTouchable.hh"

PSSDoseAveragedLET::PSSDoseAveragedLET(const G4String& name, G4int depth) : G4VPrimitiveScorer(name, depth), _fHCID(-1), _fEvtMap(nullptr)
{
    static G4ThreadLocal bool printed_ctor = false;
    if (!printed_ctor)
    {
        printed_ctor = true;
        G4cout << "[PSSDoseAveragedLET] Constructor on " << (G4Threading::IsWorkerThread() ? "WORKER" : "MASTER") << ", name=" << name << G4endl;
    }
}

PSSDoseAveragedLET::~PSSDoseAveragedLET() {}

void PSSDoseAveragedLET::Initialize(G4HCofThisEvent* HCE)
{
    _fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
    if (_fHCID < 0)
    {
        _fHCID = GetCollectionID(0);
    }
    HCE->AddHitsCollection(_fHCID, _fEvtMap);
    _fAccum.clear();

    static G4ThreadLocal bool printed_init = false;
    if (!printed_init)
    {
        printed_init = true;
        G4cout << "[PSSDoseAveragedLET] Initialize (first time this thread) on " << (G4Threading::IsWorkerThread() ? "WORKER" : "MASTER")
               << ", HCID=" << _fHCID << G4endl;
    }
}

G4bool PSSDoseAveragedLET::ProcessHits(G4Step* step, G4TouchableHistory*)
{
    const G4double edep = step->GetTotalEnergyDeposit();  // MeV
    if (edep <= 0.0) return false;

    const G4double steplen = step->GetStepLength();  // mm
    if (steplen <= 0.0) return false;

    // LET_track ~ dE/dx (MeV/mm)；数值等价于 (keV/um)
    const G4double let_track = edep / steplen;  // MeV/mm

    const G4int index = GetIndex(step);  // voxel 索引

    auto& pair = _fAccum[index];
    pair.first += let_track * edep;  // numerator   累计：Σ(LET * dE)
    pair.second += edep;             // denominator 累计：Σ dE

    return true;
}

void PSSDoseAveragedLET::EndOfEvent(G4HCofThisEvent*)
{
    // 写入 EvtMap：输出 LETd = num/den（单位：MeV/mm）
    for (const auto& it : _fAccum)
    {
        const G4int index = it.first;
        const G4double num = it.second.first;
        const G4double den = it.second.second;
        if (den > 0.0)
        {
            const G4double letd = num / den;  // MeV/mm  (== keV/um)
            _fEvtMap->add(index, letd);
        }
    }
    _fAccum.clear();
}

void PSSDoseAveragedLET::clear() { _fAccum.clear(); }

void PSSDoseAveragedLET::PrintAll()
{
    G4cout << "==> PSSDoseAveragedLET<" << GetName() << "> results (per-voxel LETd):\n";
    _fEvtMap->PrintAllHits();
}
